simulation results/check_sims.R

In BriceonWiley/IntegratedLikelihood.R:

check_sim <- function(N1, p1 = 0.1, phi1 = 0.1, theta1 = 0.1,
                           N2 = N1, p2 = 0.1, phi2 = 0.1, theta2 = 0.1,
                           alpha = 0.05, len = 1e4, odds = FALSE) {
  # browser()

  # producing and checking samples
  sample_1 <- produce_samples(N1, p1, phi1, theta1, len)
  sample_2 <- produce_samples(N2, p2, phi2, theta2, len)
  samples <- bind_cols(sample_1, sample_2) %>%
    `colnames<-`(c(
      "n100", "n101", "n110", "n111", "X1", "Y1",
      "n200", "n201", "n210", "n211", "X2", "Y2"
    ))

  # producing intervals
  Walds <-
    pmap_dfr(
      samples,
      ~ Wald_gam(
        c(..1, ..2, ..3, ..4, ..5, ..6),
        c(..7, ..8, ..8, ..10, ..11, ..12),
        alpha, odds
      )
    )
  MLEs <-
    pmap_dfr(
      samples,
      ~ {
        mle1 <- calculate_mles(c(..1, ..2, ..3, ..4, ..5, ..6))
        mle2 <- calculate_mles(c(..7, ..8, ..9, ..10, ..11, ..12))
        tibble(
          "phat1" = mle1$phat, "pihat1" = mle1$pihat,
          "lam1hat1" = mle1$lam1hat, "lam2hat1" = mle1$lam2hat,
          "phat2" = mle2$phat, "pihat2" = mle2$pihat,
          "lam1hat2" = mle2$lam1hat, "lam2hat2" = mle2$lam2hat,
          "gamhat" = log((mle1$phat * (1 - mle2$phat)) / (mle2$phat * (1 - mle1$phat)))
        )
      }
    )
  counts <-
    c(
      "N1" = N1, "M1" = N1 * 0.9, "n1" = N1 * 0.1,
      "N2" = N2, "M2" = N2 * 0.9, "n2" = N2 * 0.1, "B" = len
    )
  proportions <- c(
    "p1" = p1, "phi1" = phi1, "theta1" = theta1,
    "p2" = p2, "phi2" = phi2, "theta2" = theta2
  )
  intervals <- bind_cols(Walds)

  results <- list(
    "counts" = counts,
    "proportions" = proportions,
    "sample 1" = sample_1,
    "sample 2" = sample_2,
    "MLEs" = MLEs,
    "intervals" = intervals
  )
  results
}